namespace NBody.Cosmology
{
    using System;
    using System.Reflection;    
        
    public class Model
    {
        public event EventHandler ModelChange;
        
        private double h;
        private double omega_0;
        private double omega_l;
        private double omega_m;
        
        // Hubble constant, today, in simulation units:
        // kpc, 10^8 Msun, G = 1
        public double H0 = 4.8216e-3;

        //
        // public interface
        //
        
        public Model(double h, double omega_0, double omega_l)
        {
            this.h = h;
            this.omega_0 = omega_0;
            this.omega_l = omega_l;
            this.omega_m = omega_0 - omega_l;
            this.H0 *= h;
        }
        
        public double HParam
        {
            get { return h; }
            set { h = value; OnModelChange(); }
        }
        
        public double Omega0
        {
            get { return omega_0; }
            set { omega_0 = value; OnModelChange(); }
        }
        
        public double OmegaM
        {
            get { return omega_m; }
            set { omega_m = value; OnModelChange();}
        }
        
        public double OmegaL
        {
            get { return omega_l; }
            set { omega_l = value; OnModelChange(); }
        }
        
        public double Expansion(double z)
        {
            double t = Time(z);
            if (omega_l > 0.0)
            {
                double C = Math.Pow(OmegaM / omega_l, 1.0 / 3.0);
                double alpha = 1.5 * H0 * Math.Sqrt(omega_l);
                return C * Math.Pow(Math.Sinh(alpha * t), 2.0 / 3.0);
            }
            else
            {
                return Math.Pow(3.0 * H0 / 2.0 * t, 2.0 / 3.0);
            }
        }
        
        public double Redshift(double time)
        {
            if (omega_l > 0.0)
            {
                double c = Math.Pow(OmegaM / omega_l, 1.0/3.0);
                double alpha = 1.5 * H0 * Math.Sqrt(omega_l);
                return 1.0 / c / Math.Pow(Math.Sinh(alpha * time), 2.0/3.0) - 1.0;
            }
            else
            {
                return Math.Pow(1.5 * H0 * time, -2.0/3.0) - 1.0;
            }
        }
        
        public double Time(double z)
        {
            if (omega_l > 0.0)
            {
                double c = Math.Pow(OmegaM / omega_l, 1.0/3.0);
                double a = 1.0 / (1.0 + z);
                double alpha = 1.5 * H0 * Math.Sqrt(omega_l);
                return 1.0 / alpha * Math.Log(Math.Pow(a / c, 1.5) + Math.Sqrt(1.0 + a * a * a / c / c / c));
            }
            else
            {
                return 2.0 / 3.0 / H0 * Math.Pow(1.0 + z, -1.5);
            }
        }
        
        public double Hubble(double z)
        {
            double t = Time(z);
            if (omega_l > 0.0)
            {
                return H0 * Math.Sqrt(omega_l) / 
                       Math.Tanh(1.5 * H0 * Math.Sqrt(omega_l) * t);
            }
            else
            {
                return 2.0 / 3.0 / t;
            }
        }
            
        public double Density(double z)
        {
            double hubble = Hubble(z);
            return 3.0 * H0 * H0 / (8.0 * Math.PI) * (hubble * hubble / H0 / H0 - omega_l);
        }
        
        public double LinearGrowth(double z)
        {
            if (omega_l > 0.0)
            {
                double Om = OmegaM * Math.Pow(1.0 + z, 3.0) / (omega_l +  OmegaM * Math.Pow(1.0 + z, 3.0));
                double Ol = omega_l / (omega_l + OmegaM * Math.Pow(1.0 + z, 3.0));
                return 5.0 * Om / (2.0 * (1.0 + z)) / (Math.Pow(Om, 4.0/7.0) - Ol + (1.0 + Om / 2.0) * (1.0 + Ol / 70.0));
            }
            else
            {
                return 1.0 / (1.0 + z);
            }
        }                    
            
        public double LinearGrowthDerivative(double z)
        {
            if (omega_l > 0.0)
            {
                double Om = OmegaM * Math.Pow(1.0 + z, 3) / (omega_l +  OmegaM*Math.Pow(1.0 + z, 3));
                double Ol = 1.0 - Om;
                double A = 2.5 * Om;
                double B = 1.0 / (1.0 + z);
                double C = Math.Pow(Om, 4.0 / 7.0) - Ol + (1.0 + Om / 2.0) * (1.0 + Ol / 70.0);
                double Omdz = 3.0 * OmegaM * omega_l * Math.Pow(1.0 + z, 2.0) / Math.Pow(omega_l + OmegaM * Math.Pow(1.0 + z, 3.0), 2);
                double Oldz = -Omdz;

                // this is ddelta/dt
                double deltadot = 2.5 * Omdz * B / C - A * B * B / C - 
                                  A * B / (C * C) * ((4.0 / 7.0 * 
                                  Math.Pow(Om, -3.0 / 7.0) + 0.5 * 
                                  (1.0 + Ol / 70.0)) * Omdz + (1.0 / 70.0 * 
                                  (1.0 + Om / 2.0) - 1.0) * Oldz);
                // different from Dave, might not be correct --
                // DEFINATELY WRONG!!!!!!!!!!!!!!!!!!!!!!!!!
                double dzdt = -1.0 / Expansion(z) * Hubble(z); 
                return deltadot * dzdt;
            }
            else
            {
                //return 2.0 / 5.0 / H0 * Math.Sqrt(1.0 + z); wrong
                return H0 * Math.Sqrt(1.0 + z);
            }
        }
        
        // Standard CDM cosmology: h = 0.5, OmegaM = 1.0, BBKS transfer
        public static Model SCDM
        {
            get
            { 
                Model m = new Model(0.5, 1.0, 0.0);
                return m; 
            }
        }
        
        // current Lambda CDM -- h = 0.75 (?), OmegaL = 0.73, OmegaM = 0.26
        public static Model LambdaCDM
        {
            get
            {
                Model m = new Model(0.75, 1.0, 0.73);
                return m; 
            }
        }
        
        // from Ryden and Gunn 1987
        public static Model RG87
        {
            get
            {
                Model m = new Model(0.75, 1.0, 0.0);
                return m; 
            }
        }
        
        public static Model FromString(string cosmo_string)
        {
            Model cosmo;
            try
            {
                Type cosmo_type = typeof(Model);
                cosmo = (Model)cosmo_type.InvokeMember(cosmo_string, BindingFlags.Public | BindingFlags.Static | BindingFlags.GetProperty, null, null, null);
            }
            catch (Exception e)
            {
                Console.Error.WriteLine("Can't find cosmology model for \"{0}\"; falling back to good old SCDM", cosmo_string);
                Console.Error.WriteLine("  Problem: {0}", e.Message);
                cosmo = Model.SCDM;
            }
            
            return cosmo;
        }
        
        protected virtual void OnModelChange()
        {
            if (ModelChange != null)
                ModelChange(this, EventArgs.Empty);
        }
    }
}
